# TODO
# 1. one-sample t-test (colmean/colsd) ? but need to transpose
# 2. two-sample t-test (colmean, colvar, ncol...) ? but need to transpose
# 3. biglm or anova?
# 4. cluster t-test results???

nsplits = 12
toolpath = "/data/projects/pro42/r-nitools"
corpath = "/data/projects/pro42/cormats_4mm"

setwd(toolpath)
source("connectivity.R")
source("compare_subjects.R")

setwd(corpath)
cor.files = file.path(corpath, list.files(pattern="desc"))
cor.sublist = ni.load.cors(cor.files)

nvoxs = ncol(cor.sublist[[1]])
split.size = ceiling(nvoxs/nsplits)
split.list = seq(1,nvoxs,by=split.size)
subdist.voxlist <- foreach(i=1:nsplits, .combine="c") %dopar% 

subdist.voxlist <- ni.comparesubs(cor.sublist, start.vox=split.list[i], ninclude.voxs=split.size)


# Quick Script
library(bigmemory)
library(Rniftilib)
source("~/r-nitools/neuroimage.R")
source("~/r-nitools/base.R")

setwd("/home/data/DELAY/Nback/all_connectivity/cor_data_8mm")
list.cors = test.r2big()
list.voxmaps = test.submaps2voxmaps(list.cors)


# Functions

big.load = function(desc.file) {
    attach.big.matrix(dget(desc.file), dirname(desc.file))
}

big.rm = function(big.obj) {
    rm(big.obj)
}

r2big = function(pattern='Rda$', path='.', big.path='./bigmemory/submaps2') {
    fns = list.files(path, pattern=pattern)
    num.fns = length(fns)
    
    prefix.fns = sub(".Rda$", "", fns)
    
    list.desc = vector("character", num.fns)
    
    for (i in 1:num.fns) {
        tmp.path = file.path(big.path)
        tmp.desc = sprintf("%s/%s.desc", tmp.path, prefix.fns[i])
        
        if (file.exists(tmp.desc)) {
            cat(sprintf("Reading filename for subject #%03i out of %i\n", i, num.fns))
            list.desc[i] = tmp.desc
        } else {
            cat(sprintf("Converting subject #%03i out of %i\n", i, num.fns))
        
            load(fns[i])
            diag(cor.mat) = 1

            dir.create(tmp.path, mode="770")
            tmp.big = as.big.matrix(
                cor.mat,
                type="double",
                backingfile=sprintf("%s.bin", prefix.fns[i]),
                descriptorfile=sprintf("%s.desc", prefix.fns[i]),
                backingpath=tmp.path
            )
            rm(cor.mat)
            rm(tmp.big)
            gc()
            
            list.desc[i] = tmp.desc
        }
    }
    
    return(list.desc)
}

# Takes a list of subject descriptor files and creates a list of voxel map descriptor files
# Each descriptor file is assumed to be pointing to a file backed big matrix
# If the voxel map descriptor file does not exist, it will be created along with the big matrix
submaps2voxmaps = function(subFiles, voxPath="./bigmemory/voxmaps2") {
    
    nSubs = length(subFiles)
    
    # tmpMat should have correlations
    # with rows = # of voxels, cols = # of voxels
    tmpMat = attach.big.matrix(dget(subFiles[1]), path=dirname(subFiles[1]))
    nVoxs = nrow(tmpMat)
    rm(tmpMat)
    
    voxFiles = vector("character", nVoxs)
    
    cat("Converting", nSubs, "subject maps to", nVoxs, "voxel maps\n")
    for (v in 1:nVoxs) {
        cat(sprintf("Processing voxel #%03i...\n", v))
        
        vFilePrefix = sprintf("vox_%05i", v)    # Filename without extension
        vFile = sprintf("%s/%s.desc", voxPath, vFilePrefix)
        
        if (!file.exists(vFile)) {
            
            # Create big matrix
            # row = # of voxels, cols = # of subjects
            voxBigMat = big.matrix(
                nVoxs, nSubs,
                type = "double",
                backingfile = sprintf("%s.bin", vFilePrefix),
                descriptorfile = sprintf("%s.desc", vFilePrefix),
                backingpath = voxPath
            )
            
            # For each subject, load voxel v's vector of correlations
            for (s in 1:nSubs) {
                cat(sprintf("sub#%02i...", s))
                tmpMat = attach.big.matrix(dget(subFiles[s]), path=dirname(subFiles[s]))
                voxBigMat[,s] = tmpMat[,v]
                rm(tmpMat)
                gc()
            }
            cat("\n")
            
            rm(voxBigMat)
            gc()
        }
        
        voxFiles[v] = vFile
    }
    
    return(voxFiles)
}

# Runs kendall's W on a list of voxmaps
# i.e. looks at consistency of each voxel's connectivity map across subjects
# Each element of voxFiles should point to a big matrix descriptor file
bigKendall = function(voxFiles) {
    library(irr)

    nVoxs = length(voxFiles)
    
    # The big matrix of voxFiles[1] contains
    # rows = # of voxels, cols = # of subjects
    tmpMat = attach.big.matrix(dget(voxFiles[1]), path=dirname(voxFiles[1]))
    nSubs = ncol(tmpMat)
    rm(tmpMat); gc()
    
    cat("Computing kendall's W for", nVoxs, "voxmaps\n")
    kendallVec = vector("numeric", nVoxs)
    
    for (v in 1:nVoxs) {
        cat(sprintf("#%03i...", v))
        
        # Compute consistency of voxel maps across subjects
        voxMat = attach.big.matrix(dget(voxFiles[v]), path=dirname(voxFiles[v]))
        kendallVec[v] = kendall(voxMat[,])$value
        rm(voxMat)
        
        if ((v %% 100)==0) gc()
    }
    cat("\n")
    
    return(kendallVec)
}

# Distance-based multiple regression
# Run adonis on the between-subject comparisons of each region/voxels connectivity maps
bigAdonis = function(voxFiles, Formula, Model, ...) {
    library(vegan)
    
    Formula = sprintf("maps ~ %s", Formula)
    
    nVoxs = length(voxFiles)
    tmpMat = attach.big.matrix(dget(voxFiles[1]), path=dirname(voxFiles[1]))
    nSubs = ncol(tmpMat)
    rm(tmpMat); gc()
    
    adonisList = vector("list", nVoxs)
    
    cat("Computing across", nVoxs, "voxmaps\n")
    for (v in 1:nVoxs) {
        cat(sprintf("#%03i..", v))
        
        voxMat = attach.big.matrix(dget(voxFiles[v]), path=dirname(voxFiles[v]))
        maps = as.dist(1 - cor(voxMat[,]))
        adonisList[[v]] = adonis(as.formula(Formula), data=Model, ...)
        rm(voxMat)
        
        if ((v %% 200)==0) gc()
    }
    cat("\n")
    
    return(adonisList)
}

# Takes an adonis data structure
# and converts it into various vectors (actually matrices) of p-values
adonisGetPvals = function(adonisList) {
    
    nElems = length(adonisList)
    vecEffects = rownames(adonisList[[1]]$aov.tab)
    nEffects = length(vecEffects) - 2
    vecEffects = vecEffects[1:nEffects]
    
    pvalsMat = matrix(NA, nElems, nEffects)
    names(dim(pvalsMat)) = c("voxels", "effects")
    colnames(pvalsMat) = vecEffects
    
    for (i in 1:nElems) {
        for (j in 1:nEffects) {
            pvalsMat[i,j] = adonisList[[i]]$aov.tab$Pr[j]
        }
    }
    
    return(as.data.frame(pvalsMat))
}
